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Abstract 

We present a high order, Fourier penalty method for the Maxwell’s equations in the vicinity 
of perfect electric conductor boundary conditions. The approach relies on extending the 
smooth non-periodic domain of the equations to a periodic domain by removing the exact 
boundary conditions and introducing an analytic forcing term in the extended domain. The 
forcing, or penalty term is chosen to systematically enforce the boundary conditions to high 
order in the penalty parameter, which then allows for higher order numerical methods. We 
present an efficient numerical method for constructing the penalty term, and discretize the 
resulting equations using a Fourier spectral method. We demonstrate convergence orders of 
up to 3.5 for the one-dimensional Maxwell’s equations, and show that the numerical method 
does not suffer from dispersion (or pollution) errors. We also illustrate the approach in two 
dimensions and demonstrate convergence orders of 2.5 for transverse magnetic modes and 
1.5 for the transverse electric modes. We conclude the paper with numerous test cases in 
dimensions two and three including waves traveling in a bent waveguide, and scattering off 
of a windmill-like geometry. 

Keywords: Active penalty method. Sharp mask function, Fourier methods. Maxwell 
equations, Fourier continuation 


1. Introduction 

Pseudospectral and Fourier based methods |3E] provide a popular solution approach for 
problems involving periodic boundary conditions. Unfortunately, pseudospectral methods 
which exploit the Fourier transform do not extend easily to domains with curved bound¬ 
aries. One approach for solving partial differential equations (PDFs) on domains with curved 
boundaries is to relax the boundary condition by introducing a forcing, or penalty term to 
approximately enforce the correct boundary values. Such an approach has successfully been 
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developed for a variety of problems in fluid dynamics [SI El El E3] as well as computations 
involving turbulent flows m- Other more recent applications include using penalty equa¬ 
tions in ocean modeling [32], plasma physics [3], magneto-hydrodynamics [2E], and scalar 
advection with moving obstacles na. One signihcant drawback with such volume based 
penalty methods is the introduction of analytic errors in the penalized PDE. The resulting 
analytic error not only limits the accuracy of any numerical method, but also degrades the 
smoothness of the underlying solution. As a result of the reduced regularity in the penalized 
solution, the Fourier spectral methods typically require additional hltering steps iza. 

In recent work [33[, a new modihed penalty term was introduced to alleviate the analytic 
error due to the standard volume penalty method. The approach was examined for the heat 
and Poisson equations to obtain a third order Fourier-based method. The method was then 
extended to the Navier-Stokes equations to obtain a second order Fourier scheme. 

The focus in the current paper is on hyperbolic wave equations with an emphasis on 
Maxwell’s equations. Specihcally, we focus on the time-dependent Maxwell’s equations in free 
space in the presence of perfect electric conductors (PFC). Perfect conductors are idealized 
materials that easily conduct electricity and are accompanied with corresponding boundary 
conditions. Mathematically, PFC boundary conditions are modeled by assuming the electric 
held is normal to the boundary of the conducting material. Such a condition may then be 
converted into an appropriate Dirichlet boundary condition on the underlying PDF. 

In contrast to previous work [33] which focused primarily on elliptic and parabolic equa¬ 
tions, here a modihed approach must be applied for hyperbolic systems. Specihcally, the 
penalty term cannot be directly applied to a second order wave equation as it will introduce 
spurious oscillations in time, but rather must be introduced into the hrst order system so as 
to dampen solutions. Even with the suitable introduction of a penalty term to a hyperbolic 
system, the presence of analytic errors can signihcantly limit the accuracy of a numerical 
method. For instance we demonstrate that a conventional volume penalty method will con¬ 
verge at a rate of 0.5, namely the error scales as 0{Ax^^‘^) where Ax is the grid spacing 
of the scheme. Recent work by [7] suggests that an alternative penalization may yield hrst 
order methods, while other work BDl shows second order convergence rates for a class of hy¬ 
perbolic systems with Neumann boundary conditions. Finally, a similar in spirit approach 
[121 ESI DSl 122], where an additional penalty term is prescribed to connect subdomains, or 
to enforce boundary conditions was developed to obtain provably stable numerical schemes. 
Although that method is currently limited to low order for boundary conditions [22], it is 
hopeful that future work may lead to the development of provably higher order penalty 
methods. 

Another successful approach for solving wave problems with Fourier series is through 
Fourier extension methods [IB [22]. The methods have been very successful at obtaining 
highly accurate solutions for wave problems that do not have a divergence constraint. In 
particular, the Fourier method is combined with an iterative (alternate direction iteration) 
method to solve a sequence of elliptic problems as a means to evolve wave equations. The 
methods we propose in this paper differ as they may be discretized with an explicit in time 
method and therefore do not require solving an elliptic problem at each time iteration. 
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We emphasize that our approach is a single domain pseudospectral in space hnite dif¬ 
ference in time method. Previous single domain pseudospectral time-domain (PSTD) ap¬ 
proaches [221 E3l EH E5] cannot handle curved geometries with PEC boundary conditions. 
Our approach can be thought of as a new way to extend the single domain PSTD method 
to domains with curved geometries. In addition, our approach preserves the use of the fast 
Fourier transform (FFT) and does not suffer from dispersion errors. In subsequent devel¬ 
opments of the PSTD method US] the FFT is no longer used, multiple domains must be 
introduced, and accuracy is lost due to subdomain coupling. 

In methods such as the immersed boundary or standard penalty method, the extended 
solution is no longer smooth. The lack of smoothness then limits the convergence rate. As 
part of our approach, we ensure that the forcing creates an extension that is smooth in a 
precise sense. We demonstrate that with an appropriate modihcation and introduction of 
an active penalty term, one may achieve systematically higher order methods. Specihcally, 
we show that for problems in one dimension, one may achieve convergence rates of up to 
3.5 (the limitation currently due to time stepping), while in dimension two, one may obtain 
rates of 1.5 for transverse electric (TE) modes and 2.5 for transverse magnetic (TM) modes. 

In the hrst half of the paper we introduce the Maxwell’s equations with PEC boundary 
conditions, along with the formulation of the active penalty term. We also describe the 
analytic construction of the penalty term for TE and TM modes in dimension two. We then 
examine the analytic error in the penalty parameter for scattering of a TM mode off of a 
PEC wall. The second half of the paper focuses on the numerical implementation of solving 
the penalized Maxwell’s equations using a Fourier pseudospectral approach. Specihcally, 
we provide details on how to numerically discretize the equations in both space and time 
using equispaced grids and Fourier series. We then go on to outline details of stability 
studies in dimensions one and two and illustrate how the penalty term can be combined with 
PMLs to provide full time-dependent simulations of waves with PEC and radiating boundary 
conditions on periodic domains. In addition, we validate the approach by performing several 
numerical studies. Specihcally, we show that in dimension one, the Fourier spectral method 
does not suher from pollution (numerical dispersion) errors. We perform convergence studies 
in both one and two dimensions, showing global convergence rates of up to 3.5 in dimension 
one, 1.5 for TE modes in dimension two and 2.5 for TM modes in dimension two. Lastly, 
we illustrate the utility of the approach on some problems involving windmill shaped and 
waveguide geometries and demonstrate the natural extension to three dimensions. 

2. Basic approach 

In this paper we develop numerical Fourier methods for solving the time-dependent 
boundary value problem for Maxwell’s equations. Specihcally, we focus on solving Maxwell’s 
equations for isotropic space in the vicinity of PEC. We denote the region of isotropic space 
by Do C D = [0,D]'^, for d = 1,2,3 where [0,D]'’* is the d-dimensional cube with periodic 
boundary conditions, and the boundary P = cIDo. The Maxwell’s equations then take the 
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form 


dt 

-V X E, 

in Do X (0, T] 

(2.1a) 

dE 

~m ~ 

V X H, 

o' 

X 

O 

Cl 

(2.1b) 

V-E = 

0, 

in Do X (0, T] 

(2.1c) 

V H = 

0, 

in Do X (0, T] 

(2.1d) 

n X (E - g) = 

0, 

in F X (0, T]. 

(2.1e) 


Here we work with rescaled variables E and H so that effectively e = /i = 1 and c = 1. 
For instance, npon non-dimensionalizing Maxwell’s eqnations by rescaling {t, x) and (E, H), 
one arrives at eqnations (2.1). 

Althongh our focus will largely be on PEC, we consider a more general set of boundary 
conditions where one prescribes a Dirichlet tangential boundary condition for E as some 
general function g of time. As written in the formulation (2.1), n is the inward unit normal 
to flo, while g is the prescribed tangential component of E on F (we assume without loss of 
generality that n ■ g = 0). 

A particularly pract ical case is that of a PEC where the complimentary domain = 
[0,11]'^ \ flo (see Figure 2.1) is an electric conductor. In this case, one assumes E = H = 0 
inside fl,,. Due to the presence of either surface charges or currents, only the tangential 
component of E and normal component of H are then continuous across the interface F, 
resulting in 


n X E = 0, on F 
n • H = 0, on F. 


(2.2a) 

(2.2b) 


Note that the two boundary conditions (2.2) are equivalent. Given initial data E(x,0) = 
Eo(x), H(x, 0) = Ho(x) satisfying the compatibility conditions V ■ Eq = V • Hq = 0 and 


n X Eq = n X g, we seek a solution for (2.1) 


2.1. Penalized Equations 


We now outline how to analytically modify the equations (2.1) in the presence of PEC so 


that one may numerically solve them using time-dependent Fourier methods. The approach 
relies on extending the domain Dq to D = Dq U and suitably modifying the equations 
(2.1) inside by the introduction of a penalty term. For the practical implementation using 


Fourier methods, we take D = [0, D]'^ to be a rectangle with periodic boundary conditions. 
We then solve the full penalized equations on D with the understanding that the restriction of 
the solution to Dq represents the physical solution. Meanwhile the solution on is hctitious 
and used only to aid in the numerical computation. 
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Figure 2.1: Example domain 11 with immersed PEC with boundary F. The plot shows a 
regular grid while the points Xj, jj and yj + hiij are used in the construction of the extension 
function g. 


The modihed penalty equations take the form 


dE^ 


-V X E^, 

V X - ?7“^x/,(x)(E^ 


on n 

g), on n 
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(2.3a) 

(2.3b) 





V ■ = 0, on f2o 

V • = 0, on Qq. 


(2.4a) 

(2.4b) 


Here g is an active penalty fnnction, and Xh{'^) is a characteristic fnnction defined by 


X/x(x) 


0, for dist(x, fl^) > h 
1, for dist(x, fig) ^ h 


(2.5) 


In other words, X/i(x) = 1 if x G or within a distanc^ h to the set Qg- 

The goal is to choose g so that the penalized solntion (E^, H^) with the same initial data 
converges rapidly to the exact solntion (E, H) 


liin E^ —y E, 


lim —)• H. 


( 2 . 6 ) 


In snch a case, solving the penalized eqnations (2.3) with small h and r] provide accnrate 
approximations to the true fields. 

In the following subsections, we outline how to construct g to satisfy (2.6). Although 
there are some similarities with the procedure outlined in 


the new method described here 
differs in the sense that (at the level of a continuum PDE) the penalty term x/j(x)(E — g) is 
continuous. In other words, we choose g to continuously match E,^ at the jump discontinuity 
in Xh- In fact, when g matches m derivatives at the jump, then for a fixed p, h, and smooth 
enough boundary T, the forcing term p“^y;/i(x)(E^ — g) is C™. The construction then turns 
out to be simpler and more accurate than in 


2.2. Penalty function g in one dimension 

We start by explicitly presenting the construction of g in one dimension. For this 
construction, we assume that the boundary of the domain is located at x = 0, so that 
ffs = {x < 0} and Hq = {^^ > 0}- In addition we assume that the fields take the form 
E^ = (0, 0, E^^^(x, f)), = (0, ify^,j(x, f), 0), and that the exact solution E^(0,f) = g satis¬ 

fies a Dirichlet boundary condition at x = 0. In this case we take g = (0,0,^^(x)) to have 
only one component. The Maxwell’s equations then become 


dHy,r, 

dt 

dt 


dE^ 

dx 

dHy^r, 

dx 


V \h{x){E,^y- g,), 


(2.7a) 

(2.7b) 


where 


Xh{x) 


0, X > h 
1, X < h 


^Here dist(x,Ps) = infygna |x — y| is the distance of the point x to the set Ps. 
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The prescription is now to build gz{x) as a smooth extension of Ez^n that also goes through 
the exact boundary condition ^^(0) = g. To obtain an extension, we match m derivatives of 
Ez,n 8it X = h, for instance see Figure [2^ We choose for gz{x) to be supported on the interval 
[—L, h], where L > 0 is an order 1 parameter and h > 0 is a small parameter. Eventually 
L and h will be hxed by numerical considerations. Explicitly, we take gz{x) as a polynomial 
extension of degree (2m + 3) to smoothly extend E^^n at x = h and decay to 0 at x = —L: 

A. Matching m = 0 derivatives at x = h, which will yield a (Ax)^'^ order scheme 

gzix) = g Po,o{x) + Po,iix) (2.8) 


where 


^0,0 (a:) = 


(x — h)(x + L)‘ 


PoaM = 


x(x + LY 


hL'^ ’ ^ h{h + Ly' 

B. Matching m = 1 derivatives at x = h, which will yield a (Ax)^'^ order scheme 

gzix) = g Piflix) + Ez,r^ih) Pi,i(x) + Pi, 2 ix) (2.9) 

where 


^1,0 (a^) = 

PM = 


(x + L)^ 

~1plP~ 


[x-hf, Pi^^ix) = 


x(x + L)^ 

h(h + L)3 


1 - 


(4h + L) 
h(h + L) 


(x — h) 


x(x — h)(x + Ly 


h(h + L)3 ■ 

C. Matching m = 2 derivatives at x = h, which will yield a (Ax)^'^ order scheme 


where 


^2,0 (a;) 

^2,1 (a;) 


^2,2 (a;) 


^2,3 (a:) 


x(x + L)^ '■ 


(x + T)^(x - h)^, 
{5h + L) 


hih + iy 


1 - 


hih + L) 


(x - h) + 


x(x — h)(x + L)^ 
hih + L)^ 
x(x — h)^(x + L)‘^ 
2h(h + L)4 


(5h + L) 

hih + E) 


X 


(15^,2 + QhL + L2) 
-K) , 
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Remark 1. The important ingredient in constructing g^ is to build a smooth extension of 
Ez^r) that also satishes the exact boundary condition g. As a result, the polynomial prescrip¬ 
tion described here is not unique. In fact, other constructions - such as using an exponentially 
decaying basis [M], or solving a minimization problem - are also feasible. Future research 
involves understanding the stability properties for different extension constructions. 4|k 
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Remark 2. In practice, when using the high order extension [C] with a Fourier method, 
one only approximately computes the derivatives E'^ E'^^. Refer to details in the numerical 

implementation regarding the Fourier method. 4k 


Remark 3. The analytic convergence of the penalized solution to the underlying solution 
does not depend on the exact details of gz{x) away from the interface F. However, we 
explicitly choose gz{x) to decay to 0 at x = —L with a polynomial degree m + 2. Such a 
rate ensures that the solution E^^n is smoother at x = —L than at the point x = h. 4k 


2.3. Penalty function g for a TM mode 

In the case when the initial data Eq and Hq, and the subsequent solutions do not depend 
on the z coordinate, the components of the magnetic held decouple into a transverse magnetic 
(TM^) mode consisting of {Hx, EHy, Ef) and a transverse electric (TE^) mode consisting of 

(Ex^Ey^H^). 

In such a case, we prescribe the penalized TM^ mode to be 


OH. 


X,7] 


dt 

dH, 


yx 


dt 


dE 


z^r} 

aT 


dE^ 

dy 

dE^^y 

dx 

dHy^y 

dx 


dHx,y 

dy 


V \h{^){Ez,y- gz)- 


(2.11a) 

(2.11b) 

(2.11c) 


Here the penalty term can be taken to be g = (0,0,^^(x)), where ^^(x) depends only on 
Ez,y. We also note that since the TM^ mode only contains an E^ component, equations 
(2.11) imply that V ■ H„ = 0 and V ■ E„ = 0 for all time. 


The primary difference between the two-dimensional TM^ mode (2.11), and the one¬ 
dimensional equations (2.7) is that gz is now an extension of a two-dimensional function. To 
efficiently construct gz, we follow a similar approach to [M] where we build gz along rays 
from the boundary F. We note that in the current formulation for a fast construction of gz, 
we require that F G 

Again, we choose gz to be a continuous extension of Ez^y satisfying the exact boundary 
conditions 5^2 on x G F. To describe the construction, we make use of the following sets of 
points which are located a distance h -C 1 from F inside flo, and L away from F inside 12^: 


F -|- h = {x G IR'^ : X G Hq, dist(x, F) = h} (2-12) 

F — L = {xG]R'^:xG dist(x, F) = L}. (2-13) 


We then choose 5 '^(x) to 

(a) Match m = 0 or m = 1 normal derivatives of Ez^y at F -|- h, 

(b) Go through the exact boundary condition g{y) for any y G F, 

(c) Decay smoothly to 0 at F — L. 

The extension is constructed as follows: 












Step 1 

Build a local coordinate system surrounding the interface F in a region between T — L 
and T + h (see Figure [2^. Suppose n is the outward normal at y G F. Then one can 
write a local system (y, s) defined implicitly by 


X = y + sn(y) 


(2.14) 


where —L < s < h and y G F. 
always invert (2.14) so that s = 


For F G and sufficiently small L and h, one can 
s(x) and y = y(x) are functions of the coordinates x. 


Step 2 

Build Qzi'x.) using one-dimensional polynomials along rays. Given a point x between 
F — L and T + h, along with the corresponding point y = y(x) and distance s = s(x) 
from Step 1, the extension is 

A. Matching m = 0 derivatives at x G F -|- h, which will yield a (Ax)^'^ order scheme 


gz{^) = g{y) ^o,o(s) + hn(y)) Po,i(s) (2.15) 

B. Matching m = 1 derivatives at x G F -|- h, which will yield a (Ax)^'® order scheme 

idEz 


gz{y^) = giy) Pi,o(s) P^,r,(y + hn(y)) Pi,i(s) 


^z,ri 


. ds 


(y + ta(y)) 


P1.2M 

(2.16) 


Note that in the constructions [A] and [B], the point y-|-hn(y) is on T + h. Meanwhile, 
in construction [B], for small h 1 the expression = (n ■ y/)Ez,r] is the derivative 
of Ez^n in the normal direction n. 

Remark 4. In the case where the interface F = {x G IR'^ : ipi'x.) = 0} is described by a level 
set '0 with |V' 0 | = 1 and 


Go = {x G ; '0(x) 

> 0 } 

(2.17) 

Gs = {x G ; ipi'x.) 

< 0 }, 

(2.18) 

r + h = {x G R”^ : 0(x) 

1 = h} 

(2.19) 

F - L = {x G R”^ : 0(x) 

1 = -L}. 

( 2 . 20 ) 


In addition, n = V'0 represents a local normal to the level sets. 4^ 


Remark 5. In simple geometries, such as a circular arc, one can explicitly solve equation 
(2.14) to recover (y,s) from x. In cases where the interface F is described as the zero level 
set of a function -0 G so that '0(y) = 0 for all y G F, then one can easily recover s(x) and 
y(x) using a Newton iteration method. In this case, ip does not need to have unit normal 


(IV^I 7 ^ 1). We provide further numerical details in Section 5.2 


Remark 6. The boundary curvature cannot be inhnite. In practice, the curvature of the 
boundary F should be small enough to be resolved by the grid spacing of the spatial dis¬ 
cretization. 4 k 
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2.4- Penalty function g for a TE mode 

When the initial data and bonndary data does not depend on z, one also obtains a 
decoupled TE^ mode consisting of components Ey, Hz). Here we write the boundary 
data as g = {qx, Qy, 0), so that if n = {nx,ny, 0) is a unit normal at any point y G E on the 
boundary, the boundary condition (2.1e) reads 

!n X (E - g)) = nx{Ey - gy) - ny{Ex - gx) = 0. 


Here the penalized equations take the form 


BE. 


x,ri 


dH. 


Z,7] 


dt 

BE., 


By 


V Xh{^){Ex,r^-gx), 


y,v 


BH 


z,V 


Bt 

BH. 


Bx 


-V Xh{^){Ey^ri-gy), 


z,r) 


BE. 


XyT] 


BE. 


y,v 


Bt 


By 


Bx 


(2.21a) 

(2.21b) 

(2.21c) 


Remark 7. It is also possible to penalize only the Hz^y component of the TE^ mode using 
the equivalent Neuman PEC boundary condition (n ■ X/)Hz^y = 0. 4^ 

Remark 8. Note that in this case V ■ = 0 by virtue of the fact that only depends 

on z. Meanwhile, for a point x G Hq (or more precisely outside (P + h)), we may take the 
divergence of (2.21) to obtain 

B{V ■ E,) 


Bt 


= 0, for X outside (P + h). 


Therefore if V ■ Eq = 0, then the divergence is preserved to be zero. In the case where 
numerical spectral derivatives are used, additional care must be taken to ensure that V ■ E^ = 
0 remains zero. 4k 


The goal is to choose g to penalize the tangential component of E.,, = {Ex^y, Ey^y,0) in 
exactly the same fashion that gz penalized Ez^y in the TM^ mode. However, since there are 
now two components of g, we choose a second condition to ensure that g does not affect the 
normal component of E,, at the boundary. Namely, we choose g so that the penalty term 
penalizes only the tangential component of E^: 

n ■ (E^ - g) = 0, for y G P, 

n X (g — g) = 0 for y G P. 

The two conditions can be guaranteed provided we take g to be 


g = (E^ ■ n)n + (g - (g ■ n)n), for y G P. 

To make the construction explicit, let n = {nx,ny,0) be the normal at any point y G P 
on the boundary. Then for any point x between —L and h of P, we solve (2.14) to find 
y = y(x), s = s(x) and the corresponding normal n(y). The components of g are then 
constructed in a very similar fashion to c/z for the TM^ mode. 
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A. Matching m = 0 derivatives at x G F + h, which will yield a order scheme 


g(x) = [(E,;(y) ■ n)n + (g(y) - (g(y) ■ n)n)J Po,o{s) + E^(y + hn) Po,i{s), 
or explicitly in components 

^^(x) = [{E^^r,iy)nx + Ey^r,{y)ny)n^ + g^iy) - {gx{y)nx + gy{y)ny)nx\ Po,o(s) 

+ Ex,ri{y + hn) Po,i('S), (2.22) 

gy{x) = \^{Ex^y{y)nx + Ey^y{y)ny)ny + gy{y) - igx{y)nx + gyiy)ny)ny] Po,o(s) 

+ Ey^yiy + hn) Pqi{s). (2.23) 


It is important to note that the constructions for g^ and gy are very similar to the construction 
for gz in the TM^ mode (with the exception of having different coefficients for the Po,o term), 
and are done independently for each of the two components. Moreover, in the case of a PEC 
boundary condition, g = 0 and the penalty term g only depends on n ■ at the boundary 

r. 


2.5. Construction of g in the general case 

The more general case of constructing the extension g in higher dimensions builds on the 
general prescription described in the previous Section 2A for the TE^ mode. In particular, 
the penalty function g is chosen to penalize the tangential component of the held and to 
approximately enforce n x (E^ - g) = 0 on r. As a result, we take g to match the exact 
value of E^ at r + h, and also satisfy the tangential component of the boundary condition 

at r 


g = (E,, ■ n)n + (g - (g ■ n)n), for y G P, 
g = E,,, for y G P + h. 


Again, we may make the construction explicit. First, given any x within —L and h of P, 
solve (2.14) for y = y(x) and s = s(x). The extension is then written as 

A. Matching m = 0 derivatives at x G P + h, which will yield a (Ao:)^'^ order scheme 


gW 


(E^(y) ■ n)n + (g(y) - (g(y) ■ n)n) 


-Po,o('S) + E,,(y + hn) Pq i{s). 


Remark 9. While we have only provided the explicit construction for a (Ao:)^'^ scheme in 
the TE;j and general cases, it is possible to obtain systematically higher rates of convergence 
analytically by including additional normal derivatives in the construction of g. However, 
there are difficulties associated with obtaining stable numerical schemes in these cases. 4^ 
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Figure 2.2: Construction of g along the direction normal to the interface F. 


3. TM^ Mode Convergence Analysis 


In this section we examine the analytic convergence rate in rj and h for plane wave 
scattering solutions off of a flat wall for the two-dimensional Maxwell’s equations. The 
problem of a plane wave TM^ mode dehned on x > 0 scattering off a flat wall at x = 0 has 


a solution of the form E = {0,0, Ez{x,y,t)), H = {Hx{x,y,t),Hy{x,y,t),0). Introducing 
complex notation, Ez = "^{Ez}, = ^{Hx}, Hy = 'Sl{Hy}, the solution has components 

Ez = 


(3.1a) 

Hx = 

E^ _ gK-kx,ky)-{x,y)^ 

(3.1b) 

Hy = 


(3.1c) 


where Ei is the amplitude of the incoming wave. The PEC boundary condition .F^(0, y, f) = 0 
at X = 0 forces all of the incoming wave to be reflected back. Here we have introduced 
u ^ = kl + kl ^ 0 as the dispersion relation. We now examine the error associated with a 


solution to the penalized equations (2.11) containing the same incident wave as (3.1) with 


Xh{x,y) = 


0, X > h 
1, X < h 


(3.2) 


In this model problem, we take gz{x,y) to match the function value of at x = h with, 
for analytic simplicity, a lower order polynomial than in (2.8) given by 


x(x -1-1) 

9z{x,y) = E^^ri{h,y), 


h 


Remark 10. For simplicity, in this example we are interested in quantifying the analytic 
error induced by the improved penalty term. As a result, we take the simplest function 
gz{x,y) to be a low order polynomial which vanishes at x = —1 inside Hg. In practice, 
numerical implementations for gz{x,y) require gz to vanish more smoothly at x = —1 inside 
the obstacle region, as to avoid oscillations in the Fourier representation of <!> 
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Figure 2.3: The local construction of a ray in the vicinity of F. The function g is built as a 
polynomial along each ray. 


We now solve the penalized equations for an incoming wave with amplitude Ei and 
determine the error in the reflection due to the penalty term. For x > h, the penalized 
equations reduce to the Maxwell’s equations in free space and we can write a general solution 
as 


^z,ri 


Hx,rj 




^iu)tfi(kx,ky)-(x-h,y) _ j^^i{-k^,ky)-{x-h,y) 

* U \ 

j^ ^^iwtf^i(kx,ky)-(x-h,y) ^^i{-k^,ky)-{x-h,y) 

* 00 \ 
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(3.3) 

(3.4) 

(3.5) 



where again Ei is the amplitude of the incoming wave, while i? is a reflection coefficient to 
be determined by matching the solution across the penalty region. We note that in the exact 
unpenalized problem, Rexact = 

For X < h, we may use separation of variables and write Ez^rj = Using 

the ansatz for Ez^ni we obtain 


// = gfj( - = E^t+^kyy\ 

I o; J [zu ax J 


along with an ODE obtained from (2.11) for E(x) 
d^E 


too 


dx"^ 


+ - kl)E - —{E - g{x)) =0, x < h 


T] 


where 


^ , , X(X + 1) A/, N 

The ODE can then be simplified into the following form 


d^E 

dx'^ 


+ + A{x + = 0 


where 

and 


2 2 ; 2 -1 
7 = OJ — ky — ZOOT] 


A = ZUT] 


-1 


h{h + 1) 


E{h). 


On X < h, the ODE (3.7) has the solution 


(3.6) 


(3.7) 


E(x) = + 71 ( 27 -“ - - T'‘x^] 


(3.8) 


where 7 is the unique root with 3 fJ{* 7 } > 0 chosen to satisfy the radiation condition (expo¬ 
nential decay) for x —)■ — cxd, and Eq is the constant of integration. 

We now solve for the two unknowns Eq and R by imposing continuity of the solution Ez^y 
and Hx^n, Hy,ri at x = h for all y. Continuity of Ez,y at x = h yields 

Eo + 7l[27-^ - = y (1 - R) (3.9) 

while continuity of Hy y yields 

^yUo -f A)—y“^ — 2y“^h] = zk^Eiil + R) (3.10) 

where A = zur]~^ h(h+i) ~ equations can be used to find Eq in terms of E^ and 

R as functions of h and y. Specifically, we can eliminate Eq and write Er = REi for some 
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reflection coefficient R. Solving for R (via MATLAB’s symbolic package) and expanding in 
powers of rj and h yields 

R Ri {1 — 2ihkx) — [2\/2(l + — A)]rih + 0{h^). (3.11) 

The leading term (1 — 2ihkx) is exactly the hrst order term in the reflection coefficient 
Rexact- Therefore, R and Rexact differ by order Fixing rj = 0{Ax) and h = 0{Ax) 

yields an error of 0(Aa:^'®) in both the amplitude |-R| and phase ZR. Hence, we have a 
global error of order 1.5. 

Remark 11. Errors at both order 0(Aa;^'®) and 0{Ax‘^) appear in the expansion for |R| 
and ZR. Hence, one may initially see 2nd order convergence before observing the asymptotic 
convergence rate of 1.5. 4k 


Remark 12. One can repeat the calculation in this section by taking a static, non-active 
volume penalty term of the form Er^, where h = 0 and g = 0. Such a choice 

for a non-active penalty term r ecov ers the PEC boundary conditions, however results in a 
slow analytic convergence rate (2.6) of 0(7]^^"^). For numerical purposes, such an analytic 
convergence rate translates into a numerical scheme with global convergence 0(Ax^^^). 4k 


4. Perfectly Matched Layers (PML) 

In our current approach using Fourier methods, we work on a rectangular domain with 
periodic boundary conditions. In many applications, however, one is not interested in solving 
Maxwell’s equations in a periodic domain, but rather on an inhnite one. One major difficulty 
which arises when using a periodic computational domain to compute solutions on an inhnite 
one is the artihcial wrapping of traveling waves. Namely, waves which should radiate out on 
an inhnite domain simply wrap back into the computational domain as a result of the periodic 
boundary conditions. In this section, we outline how to eliminate the artihcial wrapping so 
that one may compute time-domain radiating solutions, such as those arising from scattering 
problems, on an ehective inhnite domain. The approach is through the introduction of a 
perfectly matched layer (PML) [D]. Although PMLs were originally introduced to eliminate 
artihcial rehections which arise from a hnite truncation of a computational domain, they are 
easily modihed to the case of a periodic domain. 

Here we outline how to modify the PML from a square domain with Dirichlet boundary 
conditions, to a periodic one. We do so for the case of a TE^ mode and note that the 
modihcation closely follows the formulation originally proposed in [5]. 

As a hrst step, we decompose the held = Hzx,q + Hzy^q into two components. In 
the absence of a PML, we choose the decomposition so that the two components evolve 
according to 

dHzx,q _ dEy q 

dt dx 

dHzy^q _ OEx^q 

dt dy 
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(4.1a) 

(4.1b) 







Although the addition of an extra equation appears redundant, the decomposition simplifies 
the resulting implementation of a PML. To add a PML we further modify the extended 
TE 2 equations (4.1) to contain an absorbing layer. For the absorbing layer, we let {axTCy) 
denote two effective material parameters. Since the domain is periodic, we simply choose the 
PML to have two bands, one vertical and one horizontal. For example. Figure 7.10 shows 
a periodic domain with two such strips outlined by dashed lines. Outside of each strip, we 
take {(Tx,(Ty) = (0,0) as a physical domain which allows for the normal propagation of the 
TF 2 mode. In such a region, one may have curved obstacles. In the PML region, we take 
o'xyO'y > 0 and modify the TF^ equations as follows: 


dEx^r, 

dt 

dEy^y 

dt 

dH^x,r) 

dt 

zy,r) 

dt 


A 

dy 

A 

dx 

A 

dx 

A 

dy 


{H^X,ri T H^yy^ 


Ey^y O'xHzXjT] 


Ex,y (^yH^yy. 


(^yEx^rj y i.EX,r] 

^xEy^y Tj x{^){Ey,ri 


9 x) 


9y) 


(4.2a) 

(4.2b) 

(4.2c) 

(4.2d) 


Through direct calculation, [S] showed that such a modification!^ results in a perfectly 
matched layer. Specifically, a wave traveling from the region where {(Jx,<Jy) = (0,0) does 
not reflect off the region where ax (or ay) is non-negative regardless of the incident angle or 
frequency. Although discontinuous jumps in ax (or ay), do not theoretically reflect waves in 
a PML, they can result in numerical reflections when computing a numerical solution. As a 
result, in practice, we choose ax (or ay) to grow linearly up to a maximum value ax,max (or 
c’'y,max)- Here, the slope and maximum value may depend on the exact problem. In practice, 
one can ramp up to the maximum value over a few wavelengths. 


5. A Numerical Fourier Algorithm 

In this section we outline the numerical method, and details we use when solving the 
penalized Maxwell’s equations. 

Let = [0, DY be the domain. Then in our scheme, we use an equispaced grid with N 
(even) points, and spacing Ax = D/N. In two dimensions we take Ay = Ax, however one 
does not in general require equal grid spacing. Grid points are denoted as 

Xj = jAx, yj = jAy, for j = 0 ,..., A - 1, 
and variables evaluated at gridpoints as uj = u{xj). 


^In the case where eo and fiQ are not 1, 
For example, see equations (2)-(3) in [3]. 


one must rescale the coefficients ax and ay in equations (4.21. 
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In the numerics, we also make use of the discrete Fourier transform of a function u{x) 
(on a domain of length D) 


N-l 


ui = T{u}= Y. 


Uj e 


—ikiXj 


(5.1) 


where 


3=0 


27r/ 

ki = -jy for 0 < Z < N/2 

ki = for iV/2 + 1 < / < iV - 1 


are the wavenumbers. The inverse is then taken as 


N-l 


U. 




^ikix 


i^j 


(5.2) 


/=0 


The discrete Fourier transform pairs also have natural extensions to higher dimensions. 


5.1. Time-stepping Details 

Linear wave equations, such as Maxwell’s equations, have an evolution governed by op¬ 
erators with purely imaginary eigenvalues. As a result, explicit time-stepping schemes may 
not be stable if the stability region does not incorporate a sufficient portion of the imaginary 
axis. The purpose of this section is to present the stability results for standard Runge- 
Kutta time stepping schemes using Fourier spectral differentiation in space in the abscence 
of penalization. 

As an example, we consider a Fourier method for the one-dimensional Maxwell’s equations 
on a periodic domain of D = 27i, given by 


dE, 

dt 

dHy 

dt 


E-Y^kE{E,}}. 


(5.3a) 

(5.3b) 


We report stability requirements for common spectral time-stepping schemes to ( |5.3 ) by 
listing the eigenvalues A to the discrete linear time evolution in Table Here the eigenvalue 
amplitude |Ap < 1 is required for stability. We denote 


r = At k EfR. 


(5.4) 


as the real parameter which combines the time step At and wavenumber k. As outlined in 
Table the simple Euler and Modified Euler schemes are always unstable. Meanwhile, RK4 
is stable provided r < y^576/72 ~ 2.83. For a d-dimensional periodic square with side length 
D, one then has k^ax = T^\/~d{N/D) with 


At < 2.83kYx 

(5.5) 

2.83 D 

(5.6) 

Tly/d N 
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Table 1: Stability for standard time-stepping schemes 


Integration of (5.3a)-(5.3b) 

Eigenvalue amplitude Ap 

Stability 

Simple Enler 

1 -)- ^2 

Unstable 

Modified Euler (RK2) 

1 + 

Unstable 

4th order Runge-Kutta (RK4) 

i 1 .^6 1 „8 

72' 57fi' 

Stable for r < 2.83 

Implicit Euler 

(l + r2)-i 

Unconditionally stable 

Integrating factor [27] 

1 

Unconditionally stable 


5.2. Solving equation (2.14) for the local coordinates 

In many applications one describes the bonndary or interface as the zero level set of 
a fnnction In this section we provide some brief nnmerical details on how one can 

nnmerically nse the level set (which may not have nnit norm) to bnild the local coordinate 


system y = y(x) and s = s(x) as the solntion to eqnation (2.14) 


We accomplish this nnmerically nsing a damped Newton method, as described in 
That is, for every grid point xj between (T — L) and (T + h), we seek the point yj on the 
zero level set of ip (i.e., the interface) snch that yj — xj is parallel to the normal direction 
n(yj) oc Vip (yj). That is, we wonld like 


/ (yt) = 


^(yj) 

(yj -Xj)x Vip (yj) 


= 0 . 


(5.7) 


Here one conld also arrive at eqnation (5.7[) by dotting and crossing (2.14) throngh with Vip 


since it is proportional to n. We note that for a Newton iteration to work, we ass nme that 
"0 G locally near the interface T so that one may compnte the Jacobian of (5.7). 


Once we have the point yj corresponding to each grid point xj, we compnte 

i’ (xj) = sigti(i/; (x,)) ||y, - x ,||2 


(5.8) 


on the grid which is now a level set fnnction with nnit norm = 1. In addition we take 
s(xj) = ip(xj) and the normal nsed in the local ray constrnction at each point is simply 
n(yi) = Vip(yj)- 

5.3. Main algorithm and details 


Discretize the spatial derivatives in (|2.3|) nsing psendospectral differentiation so that 

(5.9a) 


an. 


dt 

dEr, 

dt 


= X 

= J^-^{tkxJ^{H,}}-v-\h(x)(E,-g). 


(5.9b) 


Here, the FFT is nsed to compnte the discrete Fonrier transform (and its inverse) on 


the right hand side (RHS) of (5.9). 
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Compute g on the RHS of equation (5.9): 


(a) Build and store the local coordinate system Sj = s(xj), yj = y(xj) and normal 


n, 


= n(yj 


Do so for all grid points Xj between T — L and T + h. This may typically be done 
only once. If required, a Newton iteration with a level set may be used to solve 


equations (2.14) for each x,-. 


(b) Compute at T and T + h, i.e., at the points yj and yj + huj, respectively. 
Interpolate (via cubic interpolation) the values of at the points yj (only re¬ 
quired for the TEx or full three-dimensional cases) and yj + huj using the values 
of E„ at the equispaced gridpoints x,- 


Er,(y j) ^ Interpolate(E^) 

E??(yi + hiij) ^ Interpolate(E^) 

(c) Compute derivatives of E^ at T -|- h, i.e., at the points yj + hrvj. 

(i) For dimension one: obtain approximate derivatives E\ E" on the regular grid 
as 


B:,, = 

where c/ = 16 is a high frequency hltering parameter, followed by interpola¬ 
tion to Uj + hiij (where = ±1 in one dimension) 


^z,r]iyj + ^ Interpolate(E( ,^) 

E'z,r){yj + ^ Interpolate(E"^). 

(ii) For dimension two: obtain the required derivatives of E^, and on the 
grid 


dx 

Bn 


= 7 


= 7' 


-1 


where c/ = 16 is a high frequency hltering parameter. Interpolate the deriva¬ 
tives to the points yj + hnj G T -|- h as 


dEy 

dx 

dEr, 

By 


(jj + huj) Interpolate 
(yj -|- huj) ^ Interpolate 



(d) Using the interpolated values at yj and yj -|- hrij, use the formula [A], [B], or [C] 
from Sections 2.2-2.5 to build g(xj). 


19 













Evolve (5.9) forward in time by At using RK4 time stepping. 


Due to the spectral derivatives in (5.9), when solving either the TE^ mode or the full 
equations (2.3), a small non-zero amplitude for V ■ may arise after the 4 stages of 
RK4 (this does not occur in dimension one or for the TM^ mode when only Ez appears 
in the equations). Thus, after the 4 stages of RK4, project out the small divergence of 
E„ by computing 


V ■ E^ = • J^jEj}, (5.10) 

p(k) =(V ■ E,,)(l - X/,(x)) }, forfc=|k|^0, (5.11) 

E^ E,j — p(k)}. (5.12) 


Remark 13. Note that in part (c) of our numerical algorithm, we apply a high frequency 
hlter to obtain approximate derivatives for E^ in the construction of g. It is important to 
note that at no point do we hlter the actual solution E^ as such a procedure would destroy 
the accuracy of the algorithm. Instead, hitering E^ in the construction of g only slightly 
modihes the penalty forcing term. The parameter Cf = 16 was chosen to ensure stability of 
the numerical algorithm, while remaining small enough to preserve the overall accuracy of 
the method. ^ 

Remark 14. Note that each step of our approach makes use of well established algorithms. 
When the interface T is described by a level set, the method only requires the EFT, Newton’s 
method, bicubic interpolation, and RK4 time stepping so that implementing the method is 
straightforward given standard robust numerical packages. 


6. Stability 

In this section we discuss the stability of the numerical method in Section as well as 
the stability of the underlying penalty PDE. Specihcally, we note that there are two separate 
stability issues to consider. The hrst is the analytic stability effects that the penalty term 
has on the underlying solution, while the second is the conventional numerical stability of 
the Fourier algorithm. 


6.1. Energy and analytic stability 

In domains with PEC boundary conditions, the underlying Maxwell’s equations (2.1) 
conserve the quadratic energy 


S = l[ |E|2 + |B|2dx. 

2 Jno 

In other words, S does not depend on time. In the case of the penalized equations, the 
associated energy of the penalty held 

Sn = ]- [ |E^p + |B,,pdx 
2 JOo 
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is not exactly conserved due to the fact that no longer satishes the exact PEC boundary 
conditions. Since is close to the exact held E, the energy £r, = £ + 0{rf) for the 
appropriate 7 corresponding to the convergence rate of the method. It is important to note 
that in general £r^ could be larger (or smaller) than £, which has the interpretation of the 
penalty term pumping (or removing) a small energy into the rehected helds (see Figure 
6 . 1 ). As a result of the small increase in energy, there can be an associated weakly unstable 
eigenvalue to the penalized equations. Note that this eigenvalue can occur at the analytic 
level and is independent of numerical implementation details. Numerically the small increase 
in energy is not problematic since the errors are on the order of the numerical method. We 
also note that PMLs have similar behavior reported in the literature [III21EIDE] whose study 
is ongoing. 


Remark 15. Some numerical experiments based on varying g suggest that different formu¬ 
lations of the penalty term may act to increase or decrease the small energy difference in 
the penalized energy with respect to £. We intend to investigate the differences in future 
work. 4 k 


6.2. Stability of the numerical scheme 


Once the numerical scheme is discretized according to the algorithm in Section 5.3 


we 

examine stability by numerically computing the eigenvalues of the associated linear operators. 
One should note that the penalized equations are the sum of two operators (the wave operator 
and the penalty operator) whose eigenvalues can be independently, analytically computed. 
Unfortunately the penalty term is non-normal, so that stability is not determined by the 
eigenvalues of the penalty term alone. Alternatively, one can use energy arguments to show 
that for sufficiently small L in the penalty term g, one guarantees a strong stability preserving 
(SSP) scheme (in the Lf{VL) norm), however for such L one loses the global accuracy of the 
method. Therefore, to show stability for the current method, we compute the associated 
eigenvalues. To compute the eigenvalues, we write the numerical scheme in the form 


u 




( 6 . 1 ) 


and introduce the discretized operator 


A = 


-T] V(x)(I-G) 
-Vx 



( 6 . 2 ) 


where G is the discrete operator that approximates the penalty term with zero boundary 
condition g GE^. The Maxwell’s equations are then approximated by 


Su 

'm 


Au. 


(6.3) 


We compute the eigenvalues of A for numerous test cases and compare them to the RK4 
stability region. Specifically, we compute the eigenvalues for dimension one matching m = 
0 , 1 , 2 and for the two-dimensional TM^ mode with m = 0,1 for the domain with a hole 
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Figure 6.1: Computation of the total energy | + |B^pdx — £^(0) for a modulated 

Gaussian wave in dimension one. The penalty term forces the solution in the extended 
region, thereby adding energy to the total system. The energy in the physical domain Go, 
remains close to the initial energy in the system. After the wave reflects, the total energy is 
slightly less than the initial energy. See Section 7.1 for details associated with the specific 
scattering problem whose energy is depicted. 


removed. Although we varied different values of the parameters rj and L, we show two typical 
eigenvalue plots in Figure 6.2 indicating that the scheme is numerically stable. Finally we 
remark that when matching higher derivatives in the numerical algorithm, one needs to add 
the extra filtering step outlined in Section 5.3 part (c) for stability. Mathematically, this 
filtering step modifies the matrix G of the penalty term to make the scheme stable without 
affecting accuracy. 
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Figure 6.2: Left: Eigenvalues (dots) for the evolution operator A in one dimension (m = 0, 


N = 1024, with geometry described in Section 7.1), along with the RK4 stability region 
(line). Right: Eigenvalues for the discrete evolution operator A for the two-dimensional 


TM^ mode {m = 0, Nx = Ny = 32, with geometry described in Section 7.2). 


7. Numerical test cases 

7.1. Test 1: One-dimensional Gaussian scattering 

In this section, we perform a numerical convergence study for the active penalty method. 
We do so for the one-dimensional scattering of an incident Gaussian wave packet. Specihcally, 
we seek solutions of the form E = (0, 0, (x, f)) and H = (0, (a:, t), 0) to Maxwell’s 

equations on the domain x G [0, cx)) such that 


dH,. dE, 


dt 


dx 


dE, dH„ 


dt 


dx 


(7.1a) 

(7.1b) 
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In addition, we impose initial conditions 


Hy {x, 0 ) = Eq {f{x - xo) + f(-x - xo)} Xo(x) 
(x, 0) = Eq {f{x - Xo) - f{-x - Xo)} Xo{x) 

f{a) = e“ 2 (a) sin (cjoct) 

along with the boundary condition Ez ( 0 , t) = 0 and 


Xo (2:) 


0 X < 0 

1 X > 0 


(7.2a) 

(7.2b) 

(7.2c) 


(7.3) 


Provided xq/ct is large enough, the initial conditions simplify to a single incident Gaussian 
packet Ez (x,0) ~ Eof{x — xq). The Maxwell’s equations then have the solution 


Hy (x, t) = Eo {f{t + X - Xo) + /(t - X - Xo)} Xo{x) 
Ez (x, t) = Eo {f{t + X - Xo) - /(t - X - Xo)} Xo(a:). 


(7.4a) 

(7.4b) 


Here the hrst terms correspond to the incident wave, while the second terms correspond to 


the reflected wave. We compare the exact solution (7.4) to the numerical solution of the 
penalized equations 


dH, 


yE = :p-^{,kz,:P{Ez,y}} 


dt 

dE 


z,ri 


dt 


= 7 {ikz,7 [Hy^y]} -T] Xh ((r) {Ez,y - gz) 


with initial data (7.2), where 


Xh {x) = 


0 X > h 
1 X < h 


(7.5a) 

(7.5b) 

(7.6) 


For our test, we take Eo = I, Xo = 7, <J = l/-\/2 and Uo = 10 to be the parameters of the 
Gaussian wave packet. Meanwhile we take = [~2, 0) and the physical domain Go = [0,14) 


so that the box size is D = 16. We then integrate equations (7.5) to a final time T = 12 
using a 4th order Runge-Kutta (RK4) scheme with At = 0.2Ax, and g = h = Ax. The 
extension function 7 is constructed such that ^ 2 ( 0 ) = 0 . 

Figure [7T] compares the error for a scattered Gaussian using a non-active penalty method 
(where Qz = 0 ) to the proposed active penalty method matching m = 0 , m = 1 , and m = 2 
derivatives at the interface F = 0. For each method, we compute the asymptotic convergence 
rate and report them to be 0(Ax'^) where 7 = 0.38,1.42, 2.48 and 3.34 respectively. We note 
that the rate of approximately 1.42 when matching m = 0 derivatives is quite close to the 
predicted analytic rate of 1.5 derived in Section 

To illustrate the role of the active penalty term, we plot the penalized solution against 


the exact solution in the vicinity of the interface F. Figure 7.2 shows a standard non-active 
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Figure 7.1: Convergence rates for a scattered Guassian (7.2) with RK4 time-stepping. The 
plot shows a non-active penalty method (blue-circles), and active methods matching 0 (green- 
triangles), 1 (red-squares), and 2 (orange-diamonds) derivatives in constructing g^. The 
global convergence rates in are approximately 0.38,1.42,2.48 and 3.34 respectively. 


penalty method where = 0 for all time. Note that the poor convergence rate leads to a 
large error after the wave has reflected from the interface. Meanwhile, Figure 7.3 shows the 
penalized solution when matching m = 2 derivatives. Here, the penalty term is a smooth 
extension which matches the boundary condition at x = 0. This results in a signihcant 
increase in accuracy for the same number of grid points. 


7.1.1. One-dimensional dispersion errors 

One difficulty which arises when using finite difference methods for solving a wave equa¬ 
tion is the introduction of numerical dispersion errors. Specifically, the numerical discrete 
dispersion relation can differ from the exact analytic one at large wavenumbers. As a result, 
one must increase the resolution of the scheme, i.e., the number of grid points per wavelength 
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X X 

Figure 7.2: Comparison of the exact solution of a reflected Gaussian wave packet to that of a 
numerically computed approximation obtained using a standard non-active penalty method 
where ^2 = 0. 


(ppwl), with the wavenumber of the initial data. 

In contrast, provided one hxes the ppwl resolution, Fourier methods have been shown 
[5nj to maintain a constant error over a wide range of wavelengths. In this subsection we 
examine the pollution error for the proposed active penal ty m ethod. Here we perform the 
same test as in the previous section using the initial data (7.1) {Eq = 1, xq = 7, a = l/-\/2), 
however we vary uq G [10,500]. In the test, we fix the ppwl at either 15 or 20 so that the 
total number of grid points increases with the frequency ujq (or number of wavelengths) of 
the initial data. As in the previous test cases, we take At = CAx where C = 0.5 for m = 0 
and m = 1 derivatives and C = 0.2 for m = 2 derivatives. We also set the integration time 
T = 15, 7] = At and h = 1.002Ax. Here the factor 1.002 is taken slighly larger than 1 to 
ensure that h is at least one gridpoint away from F. In all test cases the RK4 time stepping 
scheme is used to guarantee that time discretization error is smaller than the error associated 
with the introduction of the penalty term. 

As shown in Figure 7.4, the error (in L“(r2o)) for active penalty methods remains rela¬ 


tively flat over a wide range of wavelengths. The plots also show 2nd and 4th order hnite 
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X X 

Figure 7.3: Comparison of the exact solution of a reflected Gaussian wave packet to that of 
a numerically computed approximation obtained using an active penalty method matching 
m = 2 derivatives. Note that as a result of the active penalty term Qz, the solution Ez^r] 
oscillates inside the PEC. 


difference schemes. As expected, both hnite difference schemes show an increase in error as 
the wavenumber increases. 


1.2. Test 2: Two-dimensional manufactured solutions for a domain with a circular hole 

In the following section, we test the Fourier penalty method (Section |5.3[ ) using a manu¬ 
factured solution approach on a periodic domain Q = [0,27r]^, with a circular hole removed. 
Specifically, the boundary F of the hole is given by the zero level set of the signed distance 
function 

V’ (x, y) = \l{,x- Xof + {y- yof - a. (7.7) 


The zero level set of (7.7) is circular with radius a and center {xo,yo). In our tests, we fix 
a = 2 with center xq = yo = vr. 

Here the manufactured solution approach allows for the direct convergence test of the 
penalized Maxwell’s equations. Two test problems are chosen to verify two independent 
modes of propagation supported by the two-dimensional Maxwell’s equations. We treat 
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Figure 7.4: The pollution error (in L°^(Qq)) for various active penalty methods and hnite 
difference schemes (FD2 and FD4 are 2nd and 4th order methods). For a fixed ppwl, the 
active penalty error remains flat over a wide range of wavenumbers. An active penalty 
method matching one derivative (m = 1) surpasses a 4th order finite difference scheme at 
moderate wave numbers ~ 10“^. 





























































the TM 2 mode followed by the TE;^ mode. In each case, we impose a boundary condition 
n X (E — g) = 0. Recall that when g = 0, we have the boundary condition for a PEC. 


7.2.1. TMz mode 

We seek solutions to the forced Maxwell’s equations of the form E = (0,0, {x,y,t)) 

and H = {x, y, t ), Hy (x, y, t ), 0) which satisfy 


with forcing function 
and initial conditions 


dH^ _ dE^ 

dt dy 

dHy _ dE, 

F)f 

dHy dH^ 

dt dx dy 

E = sin (x) cos (y) sin (t ), 

(x, y, 0) = 0 
Hy (x, I/, 0) = 0 
E^ (x, y, 0) = sin (x) cos (y). 


(7.8a) 

(7.8b) 

(7.8c) 


(7.9) 


(7.10a) 

(7.10b) 

(7.10c) 


The solution to (7.8) with forcing function (7.9), subject to initial conditions (7.10) is given 
by 


H^ (x, y, t) = sin (x) sin (y) sin (t) 

Hy (x, y, t) = cos (x) cos (y) sin (t) 

Ez (x, ?/, t) = sin (x) cos {y) cos (t ), 


(7.11a) 

(7.11b) 

(7.11c) 


which one can verify satishes the divergence-free criteria. 

We then compare the exact solution to the solution of the penalized equations 


dH,^y 

dt 

dHy,r, 

dt 

dEz,r, 

dt 


-E ^ {ikyE {Ez^y}] (7.12a) 

[iKE {Ez,y}} (7.12b) 

E~^ {iKE {Hy^y}] - E~^ {ikyE {Hz,^y]} - r]~^Xh {x, y) {Ez^y - Qz) + E (7.12c) 


at grid points belonging to the physical domain Hq- We note that in this instance, Qz is 
constructed to handle the non-zero boundary condition in the exact solution (7.11c). We 
integrate (7.12) to a hnal time T = l.lvr using RK4 with At = 0.4Ax, h = 2Ax, rj = 4Af, 
and L = 1. Figure 7.5| illustrates the convergence rates for the m = 0 and m = 1 cases. The 
results agree with the expected rates of 1.5 and 2.5. 
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Figure 7.5: Convergence study for the manufactured TM^ example. The global convergence 
rates in L°°{Qo) are approximately 1.5 and 2.5 for m = 0 and m = 1 respectively. Here, 
u = Ez^y), and a smoothing parameter c/ = 16 was used for the m = 1 case. No 

smoothing was required for m = 0. 


7.2.2. TEz mode 

Similarly, we may also seek solutions of the form E = (x, y, t ), Ey (x, y, t ), 0) and 

H = (0,0 , Hz (x, y, t)) over H = [0, 27r]^ such that 


dE, 

dt 

dt 

dHz 


dHz 

dy 

dHz 

dx 

dE, 

dy 


d^ 

dx 


+ E 


(7.13a) 


dt 
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(7.13b) 

(7.13c) 
















with forcing function 


and initial conditions 


E = sin (x) cos (p) sin (f), 

( 7 . 14 ) 

Ex (x, p, 0) = 0 

(7.15a) 

Ey (x, p, 0) = 0 

(7.15b) 

Hz (x, p, 0) = sin (x) cos (p). 

(7.15c) 


The solution to (7.13) with forcing function (7.14), subject to initial conditions (7.15) is 
given by 

(x, y,t) = - sin (x) sin (y) sin (t) 


Ey {x, y,t) = — cos (x) cos (y) sin (t) 


(7.16a) 
(7.16b) 

Hz{x,y,t)= sin (x) cos (?/) cos (f) (7.16c) 

which one can verify satisfies the divergence-free criteria. Using the same circular obstacle 


as in Section 7.2.1, the penalized solution is computed by integration in time (using RK4) 
of 


dE, 


x,r] 


dt 

dE, 


y,v 


dt 

OH. 


z,r) 


dt 


= {ikyE - rf-^Xh (x, y) (E^^y - Qx) 

= - V~^Xh {X, y) {Ey^y - Qy) 

= {ikyT {Ex,y}} - {ikxE {Ey^y}} + E, 


(7.17a) 

(7.17b) 

(7.17c) 


where g is constructed using g = E. That is, we penalize the electric held such that the 
tangential component at the boundary T is equal to the tangential component of the exact 


solution. Figure 7.6 illustrates the convergence rate for the m = 0 case. The parameter 


values for T, At, h, rj, and L are unchanged from Section 7.2.1 


7.3. Test 3: Solution inside a circular cavity 

We may also examine a problem similar to test case 2 where we solve the penalized 
equations on the interior of a circular cavity domain that is embedded in a periodic domain. 
In this case, we use the level set 


i’(x,y)^l - \J(x- Xaf + (y- yoY 


(7.18) 


to construct the local coordinate system for Qz. The physical, circular domain then corre¬ 
sponds to the region where 'ip > 0. 


Now consider solving (2.11) with initial conditions 


Hx,rj {p, 0, 0) = 0 
Hy,r, (p, 0, 0) = 0 
Ez,y (p, 0, 0) = Ji (a*jp) cos (i0), 


(7.19) 

(7.20) 

(7.21) 
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Figure 7.6: Convergence study for the manufactured TE^ example. The global convergence 
rate in L“(r2o) is approximately 1.5 for m = 0. Here, u = Ey^ri, Hz,n)- 


where Ji is the Bessel function of the first kind of integer order i and ajj is the positive 
real root of the order i Bessel function. The solution to the unpenalized TM^ equations is 


'I 

Hp (p, 0, t) = - Ji {aijp) sin (*0) sin [ai^t) 

CXijp 

{p, 0, t) = ^ [Ji_i (ttijp) - Ji+i (ttijp)] cos (i0) sin 
(p, 0, t) = Ji (aijp) cos (z0) cos (a*jt). 


(7.22a) 

(7.22b) 

(7.22c) 


Figure 7.7 illustrates the E^, H^, and Hy components of the solution computed using the 
proposed penalization method at time T = 0.3 with i = 6, j = 2, and agp ~ 13.5892. For 
this example, At = OAAx, p = 4At, h = 2Ax, m = 1, and L = 0.45. 


Figure 7.8 shows the convergence of the penalized solution to the exact solution (7.22). 
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Figure 7.7: Plots of the TM^ mode cavity problem solution using the penalization method. 
Note how the component tracks the penalty function Qz outside of hlo- 


7.4- Test 4-' Scattering off of a PEC cylinder 

Our last test involving the circular geometry is for a time-dependent scattering compu¬ 
tation. Specihcally, we take the following modulated Gaussian wave packet as initial data 
for a TE^ mode 


E^{x,y,0)=0 (7.23a) 

Ey{x,y,0) = ^{x - xo)e~^^^^ (7.23b) 

Hz{x, y,0) = ^{x - xo)e~^^^\ (7.23c) 

and compute the scattered wave packet off of a cylinder using our Fourier penalty method. 
Here we take the domain parameters to be G = [—1,1.5]^ with the cylinder centered at 
i^o^yo) = (0,0) with radius a = 0.2. The initial data is chosen to have a = 0.125 and 


xo = — 0 . 6 . 

To test the error, we perform the full time-dependent simulation of the scattered wave 
up to time T = 1.5080. We then compare the penalized solution with the exact analytic 
solution throughout the entire domain. The Lorentz-Mie-Debye method for electromagnetic 
scattering off of a perfectly conducting infinite cylinder is used to compute the exact solution 
for scattering of a time-harmonic plane wave in the frequency domain (see, for example, dl)- 
We then compute the time-dependent scattered solution at each grid point by taking the 
inverse Fourier transform of the exact time-harmonic solution scaled by the Fourier transform 
of the envelope of the plane wave. Since the Fourier spectrum of the Gaussian envelope is 
band-limited in finite precision, we can perform this step via inverse FFT with high accuracy. 
Figure Hi shows the convergence plot of the error, while Figure |7.10| shows a plot of the 
scattered wave. 


7.5. Test 5: Two-dimensional bent waveguide 

Next, to demonstrate some potential uses for the penalization scheme, we treat a bent 
waveguide problem. To define the waveguide geometry, we first construct five segments of a 


33 









Figure 7.8: Convergence study for the TM^ mode cylindrical PEC cavity problem. The 
global convergence rates in L°°{Qo) are approximately 1.5 and 2.5 for m = 0 and m = 1 
respectively. Here, u = and a smoothing parameter c/ = 16 was used for 

the m = 1 case. No smoothing was required for m = 0. 
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Figure 7.9: Convergence study for the TE^ mode cylindrical PEC scattering problem. The 
global convergence rate in the relative error is approximately 1.5 for m = 0. Here, u = 

^y,rii 
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(7.24) 
























Figure 7.10: Plots of the magnitudes of the (1-r) Ey^^i and components for a 

TE 2 wave scattering off of a cylinder with PEC boundary conditions. The approximate 
wavelength is A = 0.25, so that the wavelength to radius ratio is A/a = 1.25. The dashed 
lines in the plot show the region where we include a PML layer. Note that only one 
vertical and one horizontal PML strip are required. 


The constants Iq, tq, and i/os correspond to the length of straight line segments, the radii of 
circular arcs, and the y-offset for the parametrized curve respectively. We then compute a 
parametrized two-dimensional surface S (r, v) = {x (r, v) ,y {t,v) y'lp (t, v)) given by 


X 

y 


= r (r) -t- vn (r), 

-0 = C — |u| 


where c is some positive constant and 


nr = 


0 1 

-1 0 


t (t) 


with 


t(r) = 


dr 

dr 




fAr) 


(7.25) 

(7.26) 

(7.27) 

(7.28) 


The boundary of the bent waveguide corresponds to the zero level set of tp. In our example, 
we set tq = 1, Iq = 71 — 2ro, yos = 2. To avoid producing a multivalued function, we sweep 
r from 0 to 5 and v from -1 to 1 and take c = 0.5. This is sufficient for our purposes as we 
are only really interested in the signed distance function Tp in the vicinity of the zero level 
set (in particular, only signed distances h and —L away from the boundary). 

Next, we require an expression for the normal to the boundary. We hrst compute two 
tangent directions on S by differentiating with respect to r and u, and take their cross 
product. Projecting the resulting normal into the x|/-plane and normalizing to unit length 
yields 

n = sign (u) h (r) . (7.29) 
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We note that for a given r and n, we can evaluate the corresponding location {x,y), the 
level set value 'ip, and its corresponding normal n. However, in our setting, we require the 
value of Ip and its corresponding normal at a set of known locations (the grid points and a set 
of boundary points). Thus, we interpolate from equally spaced data in the rn-plane (which 
is not equally spaced in the xj/-plane) to the grid points and necessary boundary points. We 
do this once in the pre-processing stage of the algorithm before we begin our time-stepping 
scheme. 


7.5.1. TMz mode manufactured solution 

As in Section 7.2. 1 [ we hrst verify that our construction converges using a manufactured 
solution approach. In fact, we use the same manufactured solution as before and only change 
the geometry of the obstacle. In addition, due to the curvature of the waveguide boundary, 
we are required to take a smaller decay length L = 0.4 for g. Otherwise, with the exception 
of T = 0.2757r and r] = At, all other parameters are left unchanged. Figure 7.11 illustrates 
the same convergence rates as in the previous test and demonstrates the validity of the 
waveguide construction. 


7.5.2. TEz mode plane wave propagation 

More practically, consider the same waveguide geometry, but with an initial condition 
corresponding to a pulsed Gaussian. Take, for example, the initial conditions 


Ex {x,y,0) 
Ey{x,y,0) 


Hz {x,y,0) 



(7.30) 

(7.31) 

(7.32) 


with a = 0.25 and Xq = 0.5. We may then solve equations (4.2) with these initial conditions 
and the splitting Hzx{x,y,0) = Hz{x,y,0) and Hzy{x,y,0) = 0. We take ax, max = D/2Ax 
for a slab of width 0.25 and set ay = 0. In addition, T = 10, At = 0.4Ax, h = 2Ax, and 
7] = At with m = 0. In Figure 7.12 we plot a collection of snapshots of the behavior of the 
plane wave as it propagates down the waveguide. 


7.6. Test 6: Two-dimensional scattering off a windmill-like geometry 

To demonstrate scattering from objects that are not comprised of circular boundaries, we 
consider a windmill-like geometry adapted from a rhodonea curve given, in polar coordinates, 
by r = a sin {36). This trifolium is an algebraic curve corresponding to the zero level set of 

-0 (x, y) = [x^ - 1 - y^'^ — Aayx"^ -\- ay {x^ + y^'j . (7.33) 


To obtain a single smooth boundary, we shift ip and work with ip (x, y) = ip{x, y) — b with 
a = 3 and 5=1 hxed. This level set function is not a signed dist ance function (one can 
check that IV^I 7 ^ 1) so we must construct, as outlined in Section 
function ip whose zero level set coincides with that of ip. 


5.2 


a signed distance 
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Figure 7.11: Convergence study for the manufactured TM^ example applied to the bent 
waveguide geometry. The global convergence rates in L“(f2o) are approximately 1.5 and 
2.5 for m = 0 and m = 1 respectively. This agrees with our previous test case. Here, 

U (^Hx,riy Hy jjy Ez^rj') ■ 


To illustrate one possible scattering solution in the vicinity of this windmill-like geometry, 

We then 


let us consider the TM^ mode. We penalize the equations as described in Section 2.3 


add a PML to absorb outgoing scattered waves. Unlike the PML discussed in Section]^ we 
use the complex coordinate stretching interpretation of the PML [7411^135] to avoid splitting 
the penalization term in our equations. We decompose solutions to Maxwell’s equations into 
terms of the form Ez,r] {x, y) (respectively Hx,ri {x, y) and {x, y) and write 
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Figure 7.12: Six snapshots of the magnitude of the component of the plane wave prop¬ 
agating through the bent waveguide. The initial condition (t = 0) is shown in the top left 
corner, with subsequent hgures taken at times f = 2,4, 6,8,10. 
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the TM^ mode equations in the frequency domain. We then replace 


d Id 

dx 1 — iaxOJ~^ dx 
d Id 

dy 1 — iayUJ~^ dy’’ 


(7.34) 


(7.35) 


multiply both sides by the denominators, ignore all terms containing products of or ay 
with Xh (this is valid as long as the penalization and PML regions do not overlap), and 
hnally transform back to the time domain. The resulting equations are 


dH, 

dH, 


x,r] 


dE 


TT 

^ x,r) 




dt 


dE. 


z,r] 

W 

(9$ 

In 


dy 
_ dEz^y 

dx ^ 

dHx y —1 / \ / rp \ 

^x^yEz,yi 


T ^y') ^z,y T 


(7.36a) 

(7.36b) 

(7.36c) 

(7.36d) 


where <h is an auxiliary variable (initialized to zero) added to avoid integrals (terms of the 
form in the frequency domain) in the time domain representation. 

Discretization in space yields the equations 


dHx,y 

dt 

dHy,y 

dt 

dEz^y 

dt 

(9$ 

Ih 


E \lhyE t^Ez^y}} GyHx^y 

E \^ihxE \^Ez^y}} axHyy 
-E-^ {ikyE {Hx,y]] + E-^ {tkxE {Hy^y}} 

- r]~^Xh (x) {Ez,y - gz) - (cTx + o-y) Ez,y + $ 


(7.37a) 

(7.37b) 

(7.37c) 

(7.37d) 


which are integrated forward in time using RK4. We solve the problem on D = 
with initial conditions 


ffx(x,y,0)= 0 (7.38a) 

Hy(x,y,0) = -^(x-Xo)e~(^^ (7.38b) 

Ez(x,y,0)= ^(x-Xo)e~(^) (7.38c) 

which corresponds to a pulsed wave traveling in the positive x-direction. For our example, 
xo = —4 and a = 0.25. We take ax,max = o'y,max = N/2 in slabs of width 1/3. Unlike the TE^ 
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mode, we begin our simulation with both PMLs set to zero, and increase their values over 
half the duration of the simulation (we do so cubically). This is done to avoid nonphysical 
reflections of the initial plane wave off of the PML. Finally, we set T = 10, At = OAAx, 
h = 2Ax, r] = At, and L = 0.3. Table 7.13 illustrates the plane wave scattering off the 
windmill-like obstacle using m = 1 normal derivatives matched at the boundary. 

7.7. Test 1: Three-dimensional manufactured solution for a domain with a spherical hole 

To demonstrate the applicability of the method in three dimensions, we hrst test a man¬ 
ufactured standing wave solution. We seek solutions to 


(7.39a) 

(7.39b) 

(7.39c) 

(7.39d) 

(7.39e) 

(7.39f) 


OH, 

dEy 

dE, 

dt 

dz 

dy 

dHy 

dE, 

dE^ 

dt 

dx 

dz 

OH, 

dE, 

dEy 

dt 

dy 

dx 

dE^ 

dH, 

dHy 

dt 

dy 

dz 

dEy 

dH, 

dH, 

dt 

dz 

dx 

dE, 

dHy 

dH, 

dt 

dx 

dy 


with initial conditions 


E 


H 


TT 


71 


= 0 


’ 2 ^ 3 , 


= 2 (k X Eq) sin (a/S k ■ x) 


The solution to (7.39) subject to initial conditions (7.40) is given by 

E (x, t) = 2Eo cos (v^ t) cos (i/Sk • x) 

H (x,t) = 2 (k X Eq) sin {y/St) sin (A/Sk ■ x). 


(7.40a) 

(7.40b) 


(7.41a) 

(7.41b) 


which one can verify satishes the divergence-free criteria and is periodic when k = :^(1, 1,1 ) 
and Eq = (1 , —2,1). 

We solve these equations on the periodic domain fl = [0, 27r]^ with a spherical hole 
removed. Specihcally, the boundary P of the hole is given by the zero level set of the signed 
distance function 


i’ (x, y,z) = \J{x- xof + {y- yof + (z - zqY 


(7.42) 
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Figure 7.13: Six snapshots of the magnitude of the Ez^n component of a plane wave scattering 
off the windmill-like geometry. The initial condition {t = 0) is shown in the top left corner, 
with subsequent hgures taken at times t = 2,4, 6, 8,10. By the hnal snapshot, the reflected 
wave is on the verge of having completely left the computational domain. 
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with radius a and center (xq, yo, Zq). In our tests, we fix a = 2 with center xq = Po = Zq = tt. 
the penalized solution is computed by integration in time (using RK4) of 


dH^ 

dt 

dt 

OH, 

dt 

dE, 

dt 

dEy 

dt 

dE, 

dt 


{%k^E {Ey^r,}] - {tkyE {E^^y}} 

{ik^E {E^^y}} - {ik^E {E^^y}] 

{ikyT {E^,y}} - {ik^E 

[ikyE {H^^y}} - {ik^T - y~^Xh {x,y, z) {E^^y - g^) 

{ik^T {H^^y}} - E~^ {ik^E {H^^y}} - g~^Xh {x, y, z) {Ey^y - py) 

E~^ [ik^E {Hy^y}} - E~^ {ikyE {H^^y}} - g~^Xh {x, y, z) {E^^y - g^) 


(7.43a) 

(7.43b) 

(7.43c) 

(7.43d) 

(7.43e) 

(7.43f) 


where g is constructed using g = E. Figure |7.14| shows the convergence of the penalized 
solution to the exact solution at final time T = + 3 with Af = 0.4Ax, h = 2Ax, g = AAt, 

and L = 1. In three dimensions (as opposed to the TE^ mode in two dimensions), Et^, Hy, 
and Ez are not zero. As a result, we expect a slight decrease in accuracy for H as, in 
our approach, only E is penalized. This is indeed observed for the case m = 0, where the 
convergence rate for E is 1.5 and the convergence rate for H is 1. One possible improvement 
could be to design a more complicated penalization involving H to reconcile its convergence 
rate with that of E in the three-dimensional case. 


7. 8. Test 8: Three-dimensional scattering off a gyroid 

As a final example, we consider periodic scattering of a radiating dipole off of a gyroid. 


We solve the full three-dimensional Maxwell equations (7.39) on the domain O = [0,1]^ with 
initial conditions corresponding to an ideal dipole im whose initial radial envelope in the 


azimuthal plane corresponds to a Gaussian pulse of the form (7.38). The dipole is z-directed 


and lies at the point (0.2704, 0.4421, 0.3902) with II = 0.01. To compute the initial condition, 
we use the approach described in Section [7^ 

The gyroid is described by the zero level set of 

tf (x, y,z) = a — sin (27rx) cos {27iy) — sin {27iy) cos {27iz) — sin (27r2;) cos (27rx) (7.44) 


with a = 0.95 which is not a signed distance function (see Section]^). Figure 7.15 illus¬ 


trates a slice of the „ component of the radiating dipole for various times while Figure 


7.16 illustrates three level sets of the corresponding energy density | (|E^ P -|- |B,,p) of the 


computed wave solution. The figures were generated using N = 256, T = 1, At = 0.35Ax, 
g = 5Af, h = 2Ax, and L = 0.025. 
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Figure 7.14: Convergence study for the full three-dimensional Maxwell equations in the 
presence of the spherical geometry. The global convergence rate in L°°{Qo) is 1.5 for E and 
1 for H. Here, u is either or E^. 

8. Conclusions 

In this paper, we have introduced a Fourier based penalty method for solving Maxwell’s 
equations in domains with curved boundaries and perfect electric conductor boundary con¬ 
ditions. The approach relied on embedding the physical domain in a larger periodic com¬ 
putational domain, followed by the introduction of a penalty forcing term. We demonstrate 
that by constructing a penalty term that is a continuous extension of the electric held and 
that also satishes the exact boundary condition, we may systematically improve the ana¬ 
lytic convergence of the penalized PDF to the exact underlying PDF. We show by analytic 
calculations in two dimensions that one achieves high order convergence for a TM^ mode 
scattering off a straight wall. We also show through the direct computation of numerical 
eigenvalues that the scheme is numerically stable in dimension one (for m = 0,1, 2) and 
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Figure 7.15: Six snapshots of a slice of the magnitude of the component of the Gaussian 
dipole wave scattering off the gyroid. From the top left corner, the figures are taken at times 

t = 0.0999, 0.1844,0.3111, 0.4096, 0.5363,0.6348. The gyroid surface is shown in dark grey. 
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Figure 7.16: Six snapshots (matching those in Figure 7.15) of the energy density of the 
Gaussian dipole wave scattering off the gyroid. Three level sets (c = 1 in opaque blue, 
c = 0.3 in transparent cyan, and c = 0.15 in transparent yellow) of the energy density are 
shown. The gyroid surface is shown in dark ^fey. 










dimension two (for m = 0). We conclude with several numerical examples of our Fourier 
based approach. Specihcally, we show high order convergence in dimension one, as well as 
a lack of dispersion errors which typically result when solving for wave propagation at high 
frequencies. We demonstrate the approach with several more practical examples including 
propagation in a waveguide geometry and scattering off a windmill-like geometry. Finally, 
we conhrm that the method extends to three dimensions. 

Despite the simplicity of the approach, several issues can still be improved. Future work 
aims to further improve the efficiency and simplicity of constructing the extension g through 
the formulation of a minimization problem. In doing so, one can likely avoid the added step 
of solving (2.14) to compute the local coordinates. Secondly, additional stability details arise 
in dimensions two and three that are not present in dimension one and that currently limit 
the accuracy of the method to either 2.5 for TM^ modes, or 1.5 for TE^ modes. These issues 
appear due to the conditioning of the current construction for g, which relies on building 
smooth extensions along rays. The conditioning may potentially be improved by taking 
an alternative, basis based, approach to the construction of the extension g. We leave the 
investigation of alternative constructions of g for future work. Finally, one may consider a 
full (E, H) penalization. This could raise the convergence rate of H by half an order (to 
match that of E) at the expense of a more complicated scheme. 
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